#ifndef _BM_PRIMITIVE
#define _BM_PRIMITIVE
template <typename Value> class BooleanMeshVertex ; 
template <typename Value> class BooleanMeshFace ; 
template <typename Value> class BooleanMeshHalfEdge ; 
template <typename Value> class BooleanMeshEdge ; 
template <typename Value> class BooleanMesh;
template <typename Value> class IntersectionRecord ;

extern const double tiny_value ;
extern const double sq_tiny_value ;
//const double tiny_value = 1e-5 ;
const double large_value = 1e15 ;

extern std::vector<std::deque<Vector <double> > > v_curve; 
// This structure is used for merging intersection. 
extern std::map <std::pair <unsigned, unsigned>, IntersectionRecord<Value> > map_mergeInter;
	 


template <typename Value>
struct BooleanMeshPrimitive : public MeshPrimitive <Value>
{
	typedef Value            ValueType ; 
	typedef BooleanMeshFace <Value>     FaceType ;
	typedef BooleanMeshVertex <Value>   VertexType  ;
	typedef BooleanMeshHalfEdge <Value>   HalfEdgeType ; 
	typedef BooleanMeshEdge <Value>     EdgeType ;
	typedef Vector<Value>    VectorType ;
	typedef VectorType       PointType ;
	typedef BooleanMesh<BooleanMeshPrimitive<Value> > MeshType ;

};


template <typename PrimitivesType> class BooleanMesh ;

template <typename Value>
class BooleanMeshHalfEdge : public HalfEdge<Value>
{
public :

    BooleanMeshHalfEdge() 
		: HalfEdge<Value>(), mb_isallocated(false), adjacent(NULL)  {}

protected :
	typedef typename BooleanMeshPrimitive<Value>::VectorType   VectorType ;
	typedef typename BooleanMeshPrimitive<Value>::PointType    PointType ;
	typedef typename BooleanMeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef typename BooleanMeshPrimitive<Value>::MeshType     MeshType ;

public :
	bool is_allocated () { return mb_isallocated;}
	void set_allocate (bool _v = true ) {this->mb_isallocated = _v;}
	//void set_point_idx (int idx) { m_intr_point_idx = idx; } 
	//int get_point_idx () const {return m_intr_point_idx ;}
public :
	PointType m_intr_point ;
	bool mb_isallocated ;

	HalfEdgeType * adjacent ;
	MeshType *  mp_mesh; 
	friend ostream & operator << (ostream & ostr, BooleanMeshHalfEdge<Value> * he)
	{
		ostr<<he->from_v->m_coord<<' '<<he->next->from_v->m_coord<<' ';
		return ostr ;
	}

	// Array of the previous m_intr_point ;
	std::vector<PointType> m_points; 
};

template <typename Value>
struct IntersectionRecord
{
	typedef typename BooleanMeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	
	PRIMITIVE_TYPE type ;
	Segment <Value> segment ;
	Polygon <Value> polygon ;
	std::vector<std::pair<HalfEdgeType *, HalfEdgeType *> > v_merge_he; 
};

template <typename Value> 
class BooleanMeshEdge  : public Edge<Value>
{
protected:
	typedef typename BooleanMeshPrimitive<Value>::ValueType    ValueType; 
	typedef typename BooleanMeshPrimitive<Value>::FaceType     FaceType;
	typedef typename BooleanMeshPrimitive<Value>::VertexType   VertexType;
	typedef typename BooleanMeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef typename BooleanMeshPrimitive<Value>::EdgeType     EdgeType ; 
	typedef typename BooleanMeshPrimitive<Value>::VectorType   VectorType ;
	typedef typename BooleanMeshPrimitive<Value>::PointType    PointType ;

public :
	// Record the new created vertices on two meshes.
	// Build their corresponding relationship. 
	std::vector<VertexType *> v_intr_point; 

	std::vector<std::pair<PointType, VertexType *> > v_edge_point ;
private :

};
template <typename Value>
class BooleanMeshVertex : public Vertex<Value>
{
protected :
	typedef typename BooleanMeshPrimitive<Value>::ValueType    ValueType; 
	typedef typename BooleanMeshPrimitive<Value>::FaceType     FaceType;
	typedef typename BooleanMeshPrimitive<Value>::VertexType   VertexType;
	typedef typename BooleanMeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef typename BooleanMeshPrimitive<Value>::EdgeType     EdgeType ; 
	typedef typename BooleanMeshPrimitive<Value>::VectorType   VectorType ;
	typedef typename BooleanMeshPrimitive<Value>::PointType    PointType ;
	typedef typename BooleanMeshPrimitive<Value>::MeshType     MeshType ;

	typedef typename Vertex<Value>::VVIter VVIter ;
public :
    /**
	 * Calculate the area of local averaging region. 
	 */
	const ValueType get_local_area ()  const 
	{
		ValueType area = ValueType(0);
		PointType center (ValueType(0), ValueType(0), ValueType(0));
		unsigned count = 0; 


		for (VVIter iter = BooleanMeshVertex<Value>::vv_begin(); iter != BooleanMeshVertex::vv_end(); ++iter)
		{
			//Fix me later. I will define a new class Polygon. 
			center += iter->m_coord ;
			++count ;
		}
		center /= ValueType(count); 
		
		return this->m_local_area ; 
	}

	MeshType *  mp_mesh;
private :
	Value m_local_area ;

};

/**
 * This class extend the IntersectResult. 
 * It records the intersection info for constructing new faces. 
 */
template <typename Value>
struct BooleanMeshIntersectResult : public IntersectResult<Value>
{
	typedef typename BooleanMeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef BooleanMesh<Value>  MeshType; 

	// One intersection point is on one halfedge. No duplication is permitted. 
	std::vector<HalfEdgeType *> v_he; 	

	HalfEdgeType * source_he ;
	HalfEdgeType * target_he ;


	std::vector<std::vector<Vector<Value> > > cut_set ;
	std::vector<std::vector<Vector<Value> > > opposite_cut_set ;
	std::vector<std::pair<HalfEdgeType * , HalfEdgeType * > > he_set; 
	std::vector<std::pair<HalfEdgeType * , HalfEdgeType * > > opposite_he_set; 
};


/**
 * Extend the Face class
 */
template <typename Value>
class BooleanMeshFace : public Face<Value>
{

protected:
	typedef typename BooleanMeshPrimitive<Value>::ValueType    ValueType; 
	typedef typename BooleanMeshPrimitive<Value>::FaceType     FaceType;
	typedef typename BooleanMeshPrimitive<Value>::VertexType   VertexType;
	typedef typename BooleanMeshPrimitive<Value>::HalfEdgeType HalfEdgeType ; 
	typedef typename BooleanMeshPrimitive<Value>::EdgeType     EdgeType ; 
	typedef typename BooleanMeshPrimitive<Value>::VectorType   VectorType ;
	typedef typename BooleanMeshPrimitive<Value>::PointType    PointType ;
	typedef typename FaceType::FHeIter                         FHeIter; 
	typedef Face<Value>                                        ParentType ;
	typedef Polygon<Value>                                     PolygonType ;

public :
	typedef typename FaceType::FVIter                          FVIter; 
	// Array of previous points. 
	std::vector<std::pair<PointType, VertexType *> > v_face_point ;

public:

	BooleanMeshFace () : Face<Value>(), accessed (false) {}

	/**
	 * Convert the face to an oriented polygon. 
	 * @param polyogn polygon. 
	 */
	void Face2Polygon (PolygonType & _polygon)
	{
		for (FVIter iter = this->fv_begin(); iter != this->fv_end(); ++iter)
		{
			_polygon.append_vertex( (*iter)->m_coord);
		}
		
	}

	void Face2Halfedge (std::vector<HalfEdgeType *> & v_he)
	{
		v_he.clear();
		for (FHeIter iter = this->fhe_begin(); iter != this->fhe_end(); ++iter)
		{
			v_he.push_back((HalfEdgeType *) (*iter));
		}
	}

	const VectorType & operator [](unsigned _id) const
	{
		return this->m_polygon[_id] ;
	}

	
	unsigned size() const
	{
		return this->m_polygon.size();
	}
	
	/**
	 * Get the face polygon. 
	 */
	const PolygonType & get_polygon () const 
	{
		return this->m_polygon ;
	}

	/**
	 * This function will be called when constructing the mesh.
	 */
	void construct()
	{
		ParentType::construct();
		// Construct the polygon
		m_polygon.clear();
		for (FVIter iter = this->fv_begin(); iter != this->fv_end(); ++iter)
		{
			m_polygon.append_vertex( (*iter)->m_coord);
			mv_halfedge.push_back((HalfEdgeType *) (iter.get_he()));
		}
	}
	
	/**
	 *
	 */
	HalfEdgeType * get_he_by_id (unsigned idx) const 
	{
		unsigned count = 0; 		
		for (FVIter iter = this->fv_begin(); iter != this->fv_end(); ++iter, ++count)
		{
			if (count == idx) return (HalfEdgeType *)(iter.get_he());
		}
		throw ("Invalid HalfEdgeIdx");
		return (HalfEdgeType * )NULL;
	}
	/**
	 * Translate the polgyon idx to halfedge. 
	 */
	HalfEdgeType * polygon_edge_idx2halfedge(unsigned i, unsigned j) const 
	{
		HalfEdgeType * source = this->mv_halfedge[i]; 
		HalfEdgeType * target = this->mv_halfedge[j]; 

		if (source->next == target) 
		{
			return source ; 
		}
		else if (target->next == source)
		{
			return target ;
		}
		else
		{
			throw ("Invalid halfedge index.");
		}
	}


	/**
	 * Predicate the intersection of two different faces. 
	 */
	friend inline bool predicate_face_face (const FaceType * _face1, const FaceType * _face2)
	{
		return predicate_polygon_polygon(_face1->get_polygon(), _face2->get_polygon());
	}

	/**
	 * Manage the co-planar faces. 
	 */
	friend bool intersect_face_face_coplanar (const FaceType * _face1, const FaceType * _face2, BooleanMeshIntersectResult<Value> & _intr)
	{
		//std::vector<std::vector<Vector<Value> > > cut_line; 
		std::vector<pair <unsigned, unsigned > > hepair_idx ;
		std::vector<pair <unsigned, unsigned > > opposite_hepair_idx ;
		bool flag = intersect_polygon_polygon_coplanar (_face1->get_polygon(), _face2->get_polygon(), _intr, _intr.cut_set, _intr.opposite_cut_set, hepair_idx, opposite_hepair_idx);


		//std::cout<<"Cut Lines"<<std::endl;
		for (unsigned i = 0; i < _intr.cut_set.size(); ++i)
		{
			_intr.he_set.push_back(std::make_pair (_face1->get_he_by_id(hepair_idx[i].first ), _face1->get_he_by_id(hepair_idx[i].second )) );
			/*
			for (unsigned j = 0; j < _intr.cut_set[i].size(); ++j)
			{
				std::cout<<_intr.cut_set[i][j] <<' ';
			}

			std::cout<<" -- ";
			std::cout<< hepair_idx[i].first <<' '
					 << hepair_idx[i].second <<std::endl;
			std::cout<< _face1->get_he_by_id(hepair_idx[i].first ) <<' '
					 << _face1->get_he_by_id(hepair_idx[i].second ) <<std::endl;
			std::cout<<std::endl;
			*/
		}
	
		//std::cout<<"Opposite Cut Lines"<<std::endl;
		for (unsigned i = 0; i < _intr.opposite_cut_set.size(); ++i)
		{
			_intr.opposite_he_set.push_back(std::make_pair (_face2->get_he_by_id(opposite_hepair_idx[i].first ), _face2->get_he_by_id(opposite_hepair_idx[i].second ) ) );
			/*
			for (unsigned j = 0; j < _intr.opposite_cut_set[i].size(); ++j)
			{
				std::cout<<_intr.opposite_cut_set[i][j] <<' ';
			}

			std::cout<<" -- ";
			std::cout<< opposite_hepair_idx[i].first <<' '
					 << opposite_hepair_idx[i].second <<std::endl;
			std::cout<< _face2->get_he_by_id(opposite_hepair_idx[i].first ) <<' '
					 << _face2->get_he_by_id(opposite_hepair_idx[i].second ) <<std::endl;
			std::cout<<std::endl;
			*/
		}
	
		return flag; 
	}	

	friend void make_merge_key (const FaceType * _face1, const FaceType * _face2, std::pair <unsigned, unsigned> & key)
	{
		if (_face1->m_idx < _face2->m_idx) 
		{
			key.first = _face1->m_idx; key.second = _face2->m_idx; 
		}
		else
		{
			key.first = _face2->m_idx; key.second = _face1->m_idx; 
		}
	}
	/**
	 * Add the intersection result into the merge structure 
	 */
	friend void add_merge_intersection_first (const FaceType * _face1, const FaceType * _face2, BooleanMeshIntersectResult<Value> & _intr)
	{
		typename std::map <std::pair <unsigned, unsigned>, IntersectionRecord<Value> >::iterator mapIter;
		std::pair<unsigned , unsigned> key; 
		make_merge_key (_face1, _face2, key);
		IntersectionRecord <ValueType> ir ; 
		if (map_mergeInter.end() == map_mergeInter.find (key))
		{
			switch(_intr.res_type )
			{
			case SEGMENT :
				ir.type = SEGMENT ;
				ir.segment = _intr.segment ;
				// Add the halfedge 
				break; 
				;
			case POLYGON :
				ir.type = POLYGON ;
				ir.polygon = _intr.polygon ;
				break; 
				;
			default:
				throw ("Unhandled intersection result");
				;
			}
			map_mergeInter.insert(std::make_pair (key, ir));
			
		   
		}
		else
		{
		}
		//= map_mergeInter.find ();
	}
    /**
	 * Add the intersection result into the merge structure 
	 */
	friend void add_merge_intersection_second ()
	{
	}
	/**
	 * We have to rewrite the intersect part to record the face topology info. 
	 */
	friend bool intersect_face_face (const FaceType * _face1, const FaceType * _face2, BooleanMeshIntersectResult<Value> & _intr, bool & iscoplanar)
	{
		iscoplanar = false; 

		const FaceType & _pol1 = * _face1;
		const FaceType & _pol2 = * _face2;

		Vector<Value> normal1 = _pol1.normal();
		Vector<Value> normal2 = _pol2.normal();

		std::cout<<_pol1.get_polygon();
		std::cout<<_pol2.get_polygon(); 

		unsigned size1 = _pol1.size();
		unsigned size2 = _pol2.size(); 
		unsigned u0, u1, u2, u3;
		u0 =  u1 =  u2 = u3 = size1; 

		/*
		Value td0 = orient_3d (_pol1[0], _pol1[1], _pol1[2], _pol2[0]);
		Value td1 = orient_3d (_pol2[0], _pol2[1], _pol2[2], _pol1[0]);

		if (abs (td0 * td1) < tiny_value )
		{
			iscoplanar = true; 
			return intersect_face_face_coplanar(_face1, _face2, _intr);
		}
		*/

		unsigned coplanar1 = 0;

		unsigned i = 0, j = size1 - 1; 
		for (; i < size1; j = i++)
		{
			Value d0 = (_pol1[j] - _pol2[0]) * normal2 ; // Predicate
			Value d1 = (_pol1[i] - _pol2[0]) * normal2 ; // Predicate
			if (abs (d1) < Value(tiny_value)) coplanar1 ++ ;
			if (d0 * d1 <= Value(tiny_value))
			{
				if (d0 >= Value(tiny_value) && d1 < -Value(tiny_value)) { u0 = j ; u1 = i; }
				else if (d0 > Value(tiny_value) && d1 <= -Value(tiny_value)) { u0 = j ; u1 = i; }
				else if (d1 >= Value(tiny_value) && d0 < -Value(tiny_value)){ u3 = i ; u2 = j; }
				else if (d1 > Value(tiny_value) && d0 <= -Value(tiny_value)){ u3 = i ; u2 = j; }
			
			}
			
		}

		if (coplanar1 >= 3) 
		{
			//throw ("Co-planar polygons.");
			iscoplanar = true; 
			return intersect_face_face_coplanar(_face1, _face2, _intr);
		}
/*
		else if (coplanar == 2) 
		{
			throw ("Edge face intersection."); 
		}
*/
		if ( u0 == size1 || u2 == size1 ) return false;  
	
		unsigned v0, v1, v2, v3 ;
		v0 = v1 = v2 = v3 = size2; 
		unsigned coplanar2 = 0;
		for (i = 0, j = size2 - 1; i < size2; j = i++)
		{
			Value d0 = ( _pol2[j] - _pol1[0] ) * normal1 ; // Predicate
			Value d1 = ( _pol2[i] - _pol1[0] ) * normal1 ; // Predicate
			if (abs (d1) < Value(tiny_value)) coplanar2 ++ ;
			if (d0 * d1 <= Value(tiny_value))
			{
				if (d0 >= Value(tiny_value) && d1 < -Value(tiny_value)) { v0 = j ; v1 = i; }
				else if (d0 > Value(tiny_value) && d1 <= -Value(tiny_value)) { v0 = j ; v1 = i; }
				else if (d1 >= Value(tiny_value) && d0 < -Value(tiny_value)){ v3 = i ; v2 = j; }
				else if (d1 > Value(tiny_value) && d0 <= -Value(tiny_value)){ v3 = i ; v2 = j; }
			}
			/*
			if (d0 * d1 <= Value(0))
			{
				if (d0 >= Value(0) && d1 < Value(0)) { v0 = j ; v1 = i; }
				else if (d0 > Value(0) && d1 <= Value(0)) { v0 = j ; v1 = i; }
				else if (d1 >= Value(0) && d0 < Value(0)){ v3 = i ; v2 = j; }
				else if (d1 > Value(0) && d0 <= Value(0)){ v3 = i ; v2 = j; }
			}
			*/
		}
		if (v0 == size2 || v2 == size2 ) return false; 

		Value d0 = orient_3d (_pol1[u0], _pol1[u1], _pol2[v0], _pol2[v1]);
		Value d1 = orient_3d (_pol1[u2], _pol1[u3], _pol2[v3], _pol2[v2]);
	
		if (!(d0 <= 0 && d1 <= 0)) return false;

		Value d2 = orient_3d (_pol1[u0], _pol1[u1], _pol2[v3], _pol2[v2]);
		Value d3 = orient_3d (_pol1[u2], _pol1[u3], _pol2[v0], _pol2[v1]);

		// Construct the intersection here. 
		std::cout<<"coplanar "<<coplanar1<<' '<<coplanar2<<std::endl;
		Vector<Value> source, target; 
		unsigned u, v; 
		HalfEdgeType * he; 

		// The intersection part is a point. 
		//if (d0 == 0 || d1 == 0)
		if ((d0 > -tiny_value && d0 < tiny_value ) || (d1 > -tiny_value && d1 < tiny_value ))
		{
			std::cout<<"POINT Degenerated cases."<<std::endl;
			//throw ("POINT Degenerated cases. ");
			//std::cout << _intr.v_he[_intr.v_he.size() - 1]->face->m_idx << std::endl;
			if (d2 < Value(0))
			{
				if (u1 == u2) u = u3; else u = u2; 
				source = intersect_seg_triangle (_pol2[v3], _pol2[v2], _pol1[u0], _pol1[u1], _pol1[u]); 
			}
			else 
			{
				if (v1 == v2) v = v3; else v = v2;
				source = intersect_seg_triangle (_pol1[u0], _pol1[u1], _pol2[v0], _pol2[v1], _pol2[v]);
			}
			_intr.res_type = POINT ;
			_intr.point = source ;
			return true; 
		}


		if (d2 < -Value(tiny_value))
		{
			if (u1 == u2) u = u3; else u = u2; 
			source = intersect_seg_triangle (_pol2[v3], _pol2[v2], _pol1[u0], _pol1[u1], _pol1[u]); 
			he = _face2->polygon_edge_idx2halfedge(v3, v2); 
			_intr.source_he = he; 
		}
		else 
		{
			if (v1 == v2) v = v3; else v = v2;
			source = intersect_seg_triangle (_pol1[u0], _pol1[u1], _pol2[v0], _pol2[v1], _pol2[v]);
			he = _face1->polygon_edge_idx2halfedge(u0, u1);
			_intr.source_he = he; 
		}


		if (d3 < -Value(tiny_value))
		{
			if (u1 == u2) u = u3;  else u = u2; 
			target = intersect_seg_triangle (_pol2[v0], _pol2[v1], _pol1[u0], _pol1[u1], _pol1[u]);
			he = _face2->polygon_edge_idx2halfedge(v0, v1);
			_intr.target_he = he; 
		}
		else
		{
			if (v1 == v2) v = v3; else v = v2; 
			target = intersect_seg_triangle (_pol1[u2], _pol1[u3], _pol2[v0], _pol2[v1], _pol2[v]);
			he = _face1->polygon_edge_idx2halfedge(u2, u3);
			_intr.target_he = he; 
		}

		if ((source - target).square_module()< sq_tiny_value) 
		{
			_intr.res_type = POINT ;
			_intr.point = source; 
			return true; 
			//throw (" Point . Degeneracy cases. ");
		}

		// The intersection part is a segment. 
		_intr.res_type = SEGMENT ;
		_intr.segment = Segment<Value>(source, target);

		/*
		if (coplanar1 != 2 && coplanar2 == 2)
		{
			return true ;
		}
		else return false; 
		*/
		return true; 
	}


	/**
	 * This function will be called when updating data of mesh. 
	 */
	void update ()
	{
		ParentType::update();
		m_polygon.clear();
		for (FVIter iter = this->fv_begin(); iter != this->fv_end(); ++iter)
		{
			m_polygon.append_vertex( (*iter)->m_coord);
		}
	}
protected :
	PolygonType m_polygon ;
	std::vector<HalfEdgeType *> mv_halfedge;  

public :
	bool accessed ;

public :
//	bool splitted ; // True : has been splitted, so cannot be copied; False : otherwise; 
};


#endif
